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Abstract. We examine the application of a new implicit unconditionally- 
stable high-resolution TVD scheme to steady-state calculations. It is a 
member of a one-parameter family of explicit and implicit second-order ac- 
curate schemes developed by Harten for the computation of weak solutions 
of one-dimensional hyperbolic conservation laws. This scheme is guaranteed 
not to generate spurious oscillations for a nonlinear scalar equation and a 
constant coefficient system. Numerical experiments show that this scheme 
not only has a fairly rapid convergence rate, but also generates a highly- 
resolved approximation to the steady-state solution. A detailed implemen- 
tation of the implicit scheme for the one- and two-dimensional compressible 
inviscid equations of gas dynamics is presented. Some numerical computar 
tions of one- and two-dimensional fluid flows containing shocks demonstrate 
the efficiency and accuracy of this new scheme. 


§1. Introduction 

Conventional schemes for the solution of nonlinear hyperbolic conservation 
laws are linear and L 2 -stable when considered in the constant coefficient case 
[1]. There are three major difficulties in using such schemes to compute 
discontinuous solutions of a nonlinear system, such as the compressible Euler 
equations: 

(i) Schemes that are second (or higher) order accurate may produce oscil- 
lations wherever the solution is not smooth. 

(ii) Nonlinear instabilities may develop in spite of the /^-stability in the 
constant coefficient case. 


fResearch Scientist, Computational Fluid Dynamics Branch. 
^Associate Professor, School of Mathematical Sciences. 



(iii) The scheme may select a nonphysical solution. 

It is well known that monotone conservative difference schemes always 
converge and that their limit is the physical weak solution satisfying an 
entropy inequality. Thus monotone schemes are guaranteed to not have 
difficulties (ii) and (iii) above. However, monotone schemes are only first-order 
accurate. Consequently, they produce rather crude approximations whenever 
the solution varies strongly in space or time. 

When using a second (or higher) order accurate scheme, some of these 
difficulties can be overcome by adding a hefty amount of numerical dissipation 
to the scheme. Unfortunately this process brings about an irretrievable 
loss of information that exhibits itself in degraded accuracy and smeared 
discontinuities. Thus, a typical complaint about conventional schemes which 
are developed under the guidelines of linear theory is that they are not robust 
and/or not accurate enough. 

To overcome the above difficulties, we consider a new class of schemes that 
is more appropriate for the computation of weak solutions (i.e., solutions 
with shocks and contact discontinuties) of nonlinear hyperbolic conservation 
laws. These schemes are required (a) to be total variation diminishing in the 
nonlinear scalar case and the constant coefficient system case (see[2,3]), and 
(b) to be consistent with the conservation law and an entropy inequality [4-6]. 
The first property guarantees that the scheme does not generate spurious 
oscillations. We refer to schemes with this property as Total Variation 
diminishing ( TVD ) schemes (or total variation nomncreasing - TVNI [2]). 
The latter property guarantees that the weak solutions are physical ones. 
The motivation for considering this class of schemes lies in the convergence 
theorem that states: “If a finite difference scheme is consistent with the 
conservation law and its entropy inequality, and the scheme is stable in the 
sense of uniformly bounded total variation in x, then the scheme is convergent 
and its limit is the physical weak solution of the initial value problem” (see 
[3]). As a result of this theorem, schemes in this class are guaranteed to avoid 
difficulties (i) to (iii) mentioned above. 

The class of TVD schemes contains that of monotone schemes, but 
is significantly larger, as it includes also second-order accurate schemes. 
Existence of second-order accurate TVD schemes was demonstrated in [2,3,7- 
9]. Unlike monotone schemes, TVD schemes are not automatically consis- 
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tent with the entropy inequality. Consequently, some mechanism may have 
to be explicitly added to a TVD scheme in order to enforce the selection 
of the physical solution. In [2] and [10], Harten and Harten and Hyman 
demonstrate a way of modifying a TVD scheme to be consistent with an 
entropy inequality. 

Second-order accurate TVD schemes are genuinely nonlinear in the 
sense that the schemes are necessarily nonlinear even in the constant 
coefficient case. Consequently, nonstandard design principles are needed 
for the construction of second-order accurate TVD schemes. In [2], Harten 
develops a rather general technique for converting a first-order accurate ex- 
plicit TVD scheme into a second-order accurate one. The second-order ac- 
curacy is obtained by applying the first-order accurate TVD scheme to an 
appropriately modified flux; the so modified scheme remains TVD. This is 
due to the fact that the TVD property is independent of the particular form 
of the flux as long as it has bounded modified characteristic speeds. 

In [11,12], we have examined the application of an explicit second-order 
accurate TVD scheme [2] to steady-state calculations. Numerical experi- 
ments show that this explicit scheme generates nonoscillatory, highly accurate 
steady-state solutions. 

The time-consistent approach to steady-state may be thought of as the 
temporal decay of the initial error. The rate of decay associated with the 
pure initial value problem is 0(t~ 1 / 2 ) in the infinite domain and 0(t~ *) in 
the periodic case (see [13,14]). Therefore the use of an explicit scheme in a 
time-consistent approach may turn out to be extremely expensive. 

To retain the characteristic of highly-resolved steady-state solutions by 
explicit second-order accurate TVD schemes without the disadvantage of 
slow convergence rate of explicit schemes, we considered in [11] the following 
two possibilities: (1) First obtain an approximation to the steady-state by 
using a conventional implicit scheme, and then use the second-order accurate 
TVD scheme as a “post-processor". (2) Use a first-order accurate implicit 
scheme in delta-formulation and replace the explicit operator by the explicit 
second-order accurate TVD scheme. 

We have found (in one dimension) that both these strategies reduce the 
overall computational effort needed to obtain the steady-state solution of the 
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explicit second-order accurate TVD scheme. Alternative (1) is a possible way 
of speeding up the convergence process by providing a better initial condition 
for the explicit second-order accurate TVD scheme. Alternative (2) can be 
viewed as a relaxation procedure to the steady-state solution. Numerical 
experiments of [11] show that the computational effort is not drastically 
decreased, although the stability limit is higher than the explicit counter- 
part. 

Recently, Harten [3] has extended the class of explicit TVD schemes to 
a more general category which includes a one-parameter family of implicit 
second-order accurate schemes. Included in this class are the commonly used 
time-differencing schemes such as the backward Euler and the trapezoidal 
formula. 

This paper is a sequel to [11]. In here we investigate the application 
to steady-state calculations of this newly developed implicit second-order 
accurate scheme that is unconditionally TVD. This scheme is guaranteed not 
to generate spurious oscillations for nonlinear scalar equations and constant 
coefficient systems. Numerical experiments show that this scheme has a rapid 
convergence rate, in addition to generating a highly-resolved approximation 
to the steady-state solution. We remark that all of the analysis on the new 
scheme is for the initial value problem. The numerical boundary conditions 
are not included. 

In the present paper we stress applications rather than theory, and we 
refer the interested reader to [2,3] for more theoretical details. In the next 
section, we will briefly review the notion of TVD schemes and describe the 
construction of the second-order accurate TVD scheme from a first-order 
accurate one for scalar one-dimensional hyperbolic conservation laws. The 
generalization to one-dimensional hyperbolic systems will be described in 
section IH. A description of the algorithm and numerical results for the one- 
and two-dimensional compressible inviscid equations of gas dynamics will be 
presented in sections IV and V. 

The paper is organized as follows: 


I. Introduction 

H. TVD Schemes for One-Dimensional Nonlinear Scalar Hyperbolic 
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Conservation Laws 

2.1 Explicit TVD schemes 

2.2 Implicit TVD Schemes , 

2.3 First-Order Accurate Backward Euler Implicit TVD Scheme 

2.4 Conversion to Second-Order Accurate Schemes 

2.5 Enhancement of Resolution by Artificial Compression 

2.6 Linearized Version of the Implicit TVD Scheme 

HI. Generalization to One-Dimensional Hyperbolic Systems of Conservation 
Laws 

IV. Applications to One-Dimensional Compressible Inviscid Equations of 
Gas Dynamics 

4.1 Description of the Implicit Algorithm 

4.2 Numerical Example 

4.3 Discussion of Numerical Results 

V. Applications to Two-Dimensional Compressible Inviscid Equations of 
Gas Dynamics 

5.1 N um erical Fluxes in Two-Dimensions 

5.2 Extension of the Explicit TVD Scheme by the Fractional Step 

Method 

5.3 Extension of the Implicit TVD Scheme by the Alternating Direction 

Implicit (ADI) Method 

5.4 Numerical Example 

5.5 Discussion of Numerical Results 

VI. Concluding Remarks 


§2. TVD Schemes for One-Dimensional Scalar Hyperbolic Conservation Laws 

Several techniques for the construction of nonlinear, explicit, second-order 
accurate, high- resolution, entropy satisfying schemes for hyperbolic conser- 
vation laws have been developed in recent years. See, for example, van Leer 
[7], Colella and Woodward [15], Harten [2], and Osher and Chakravarthy 
[9]. We can also view these schemes as shock-capturing algorithms based 
on either an exact or approximate Riemann solver. From the standpoint of 
numerical analysis, these schemes are TVD for nonlinear scalar hyperbolic 
conservation laws and for constant coefficient hyperbolic systems. Entropy 
satisfying TVD schemes have the property that they do not generate spurious 
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oscillations and that the weak solutions are physical ones. The goal of con- 
structing these highly nonlinear schemes is to simulate complex flow fields 
more accurately. TVD schemes are usually rather complicated to use com- 
pared with the conventional shock-capturing methods such as variants of the 
Lax-Wendroff scheme. 

In [3], Harten introduced the notation of implicit TVD schemes. To keep 
this paper somewhat self-contained, we will review the construction of the 
backward Euler TVD schemes for the initial value problem. This is the only 
unconditionally stable TVD scheme belonging to the one-parameter family 
of TVD schemes considered in [3]. Although the trapezoidal formula belongs 
to the same one-parameter family, the CFL like requirement for this implicit 
scheme to be TVD is 2; thus it is neither efficient nor practical for steady- 
state calculation. Before we proceed with the description of the construction, 
we will first give preliminaries on the definition of explicit and implicit TVD 
schemes and show a few examples. 


§2.1 Explicit TVD Schemes 


Consider the scalar hyperbolic conservation law 


du df(u) 
dt + dx 



( 2 . 1 ) 


where a(u) — df/du is the characteristic speed. A general three-point 
explicit difference scheme in conservation form can be written as 

u y"*" 1 = M y ~ M/j-t-i/2 — f j— 1 / 2)1 (2-2) 

where f j+i / 2 = /(«”, u ]+i)’ ^ = At/ Ax, with At the time step, and Ax the 
mesh size. Here, = u(jAx,nAt) and / is a numerical flux function. We 

require the numerical flux function / to be consistent with the conservation 
law in the following sense 


f(Uj, Uj ) = f(Uj) 


(2.3) 
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Consider a numerical scheme with numerical flux functions of the following 
form 


! 1+1/2 ~ M + fj+l — Q( a i+1/2)&1+1/2U] 


(2.4) 


where fj = f(uj), Aj+i/ 2 u = Uj+i — Uj and 


r 

a j+ 1/2 = ' 


(fj- 1-1 “ fj)/ A j+ 1/2 U 
a(uj ) 


Aj-|-1/2^ ^ 0 
&j+l/2 u — 0. 


(2.5) 


Here Q is a function of aj+i /2 and X. The function Q is sometimes referred 
to as the coefficient of numerical viscosity. Figure (2.1) shows some examples 
for the possible choice of Q. Three familiar schemes with the numerical fluxes 
of the form (2.4) are: 

(a) A form of the Lax-Wendroff (L-W) scheme with 

fj+ 1/2 — 2 \.fj + /r-H “ M a >+i/ 2 ) 2 Ay+i/ 2 tt] (2.6) 


where Q{a j+l/2 ) = X(a J+ i /2 ) 2 

(b) Lax-Friedrichs (L-F) scheme with 

fj+1/2 = jjrl/j'H" ■^’+ 1 — ^>+1/ 2 W ] (2-7) 

where Q(aj+i/ 2 ) — 1 

(c) A generalization of the Courant-Isaacson-Rees (GCIR) scheme with 

7 3+1/2 = ^[fj + fj+1 — K+i/2|A j+ i/ 2 u] ( 2.8) 

where Q(°j- {- 1/ 2 ) = |®y-j-i/2| 

We define the total variation of a mesh function u to be 
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(2.9) 


oo oo 

TV(u) — £ K+i - u j\ = £ l A i+i/2«l 

j— — oo j— — oo 

We say that the numerical scheme (2.2) is TVD if 

TV(u n+l ) < TV{u n ) (2.10) 

It can be shown that a sufficient condition for (2.2), together with (2.4), to 
be a TVD scheme is [2] 


j+l/2 ~ ^[— a J+l/2 + Q(Oj + 1/2)] > 0 

(2.11a) 

xc f+ 1/2 = g [ fl J+i/2 + Q{aj+ 1 / 2 )] > 0 

(2.11b) 

M^y- f. 1/2 ~h = 1 / 2 ) < 1 

(2.11c) 


Applying condition (2.10) or (2.11) to the above three examples, it can be 
easily shown that the L-W scheme is not a TVD scheme, and the latter two 
schemes are TVD schemes. Note that there is a further distinction between 
the L-F scheme and GCDR scheme: the L-F scheme is consistent with an 
entropy inequality whereas the GCIR is not [6]. 

It should be emphasized that condition (2.11) is only a sufficient condition; 
i.e., schemes that fail this test might be still TVD. The L-W scheme, besides 
failing condition (2.11), does not satisfy (2.10). 


§2.2 Implicit TVD Schemes 

Now we consider a one parameter family of three-point conservative 
schemes of the form 

u j'^ 1 + 2 — f j—\j 2 ) — u1 j ~ Ml — v)(f j-\- 1/2 f j— 1 / 2 ) (2-12) 


7 



where rj is a parameter, X = At/ Ax, fj+1/2 = /( w i> w i+i)> fj+1/2 ~ 
and is the numerical flux (2.4). This one- 

parameter family of schemes contain implicit as well as explicit schemes. 
When T) = 0, (2.12) reduces to (2.2), the explicit method. When rj ^ 0, 
(2.12) is an implicit scheme. For example: if r/ = 1/2, the time differencing is 
the trapezoidal formula and if rj = 1, the time differencing is the backward 
Euler method. To simplify the notation, we will rewrite (2.12) as 


L-u n+ 1 =R-u n (2.13) 

where L and R are the following finite-difference operators: 

(L • u) j = Uj + Xt/(7 j+ i / 2 ~ 1 3 - 1 / 2 ) (2.14a) 

(R • u) j = Uj - X(1 - v)Uj+i /2 ~ ~f j — 1/2) (2.14b) 

A sufficient condition for (2.12) to be a TVD scheme is that 

TV(R • v ) < TV{v ) (2.15a) 

and 

TV(L ■ v) > TV{v) (2.15b) 

A sufficient condition for (2.15) is the CFL-like condition 

|^ctj+i/2| < ^Q(Oj+ 1/2) < " _ ■ " (2.16) 


where aj+1/2 is defined in equation (2.5). For a detailed proof of (2.15) and 
(2.16), see [3]. Observe that the backward Euler implicit scheme, 77 = 1 
in (2.12) is unconditionally TVD, while the trapezoidal formula, i) = 1/2 
is TVD under the CFL-like restriction of 2. The forward Euler explicit 
scheme, r\ — 0 or (2.2), is TVD under the CFL condition of 1. We remark 
that three-point conservative TVD schemes of the form (2.12) are generally 
first-order accurate in space. When rj — 1/2, it is second-order accurate in 
time. 
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§2.3 First-Order Accurate Backward Euler Implicit TVD Scheme 

In this paper, we are only interested in efficient high resolution time depen- 
dent methods for steady-state calculations. The backward Euler implicit 
TVD scheme is the best choice in this one-parameter family of TVD schemes. 
In this section we will review the proof that the backward Euler scheme is 
unconditionally TVD. In section 2.4, we will decribe the technique of con- 
verting the first-order accurate unconditionally TVD scheme (2.12) with rj = 
1 into a second-order accurate one. 

The backward Euler three-point scheme in conservative form can be written 
as 


*" +l + Mfjti/2 - 77—1/2) = «" (2.17a) 


with 


f{Uj>Uj+ 1) — 2 i/j /? + ! Q( a j-H/2)Aj-|-l/2w] 


(2.17b) 


where 


a i+ 1/2 = \ 


if 3 + 1 fj)/ A/'+ 1/2^ 


A7+1/2M ^ 0 
A j+ i/ 2 « =0. 


(2.17c) 


For the purpose of this paper, the function Q(z) in (2.17b) is chosen to be 



£(z 2 /5 + 8) 


M < 6 

\z\ > 8 . 


(2.17d) 


which is a nonvanishing, continuously differentiable approximation to \z\. 
This choice of Q is the least dissipative among the class of three- point entropy 
satisfying TVD schemes that we are considering (see figure (2.1)). Note that 
if 6 = 0 in (2.17d), (i.e., Q(z) = |z|) then (2.17a) is a first-order accurate, 
upstream differencing, backward-Euler scheme. In order to see the above 
(with a more familiar notation for upstream differencing representation) and 
to guide us for later development, we introduce the following notation 

<?%)= i[«(2)±*] (2.18a) 
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and note that 


CHz) >0 for all 2 (2.18b) 

Using (2.17c) and (2.18a), we can rewrite the numerical fluxes fj ± 1/2 in 
(2.17b) as 

fj+ 1/2 = fj ~ C'~(aj+i/ 2 )^j+i/ 2 W (2.19a) 

f j— 1/2 = fj i/ 2 )Aj_ 1 / 2 U (2.19b) 

It follows from (2.19) that (2.17) can be written in the form 

«? +l - XC-(ay+ 1 1 /2 )A, +1/2 «"+ 1 + XC+(«J>+ 1 1 /2 )A,-_i /2 ii’ , + 1 = u] 

( 2 . 20 ) 

Now, if 8 — 0 in (2.17d), then C^(z) — ^{\z\/tz), and (2.20) is a first-order 
accurate, upstream differencing, backward Euler implicit scheme. Equation 
(2.17) differs from the upstream spatial differencing (with 8 — 0) by the 
addition of a so called numerical viscosity term with a coefficient 6 > 0; 8 has 
the dimensions of velocity. The addition of this viscosity term accomplishes 
two things: (i) it makes the numerical flux continuously differentiable, (ii) it 
rejects stationary expansion shocks (see [2] for more details). In general, a 
larger 8 brings about improvement in entropy enforcement at the expense of 
some loss of resolution. However, from numerical experiments the dependence 
on 6 is rather slight with 6 = 0.025,0.05 and 0.125. 

We show now that C±(z) > 0 implies that the scheme (2.17) is uncondi- 
tionally TVD (i.e., condition (2.10) is satisfied), independent of the value of 
X = At/ Ax in (2.17a). 

To see that, we subtract (2.20) at j from (2.20) at j -j- 1 and get after 
rearranging terms that 


1 + XC- +1/2 + XC+ +1/2 jA, + 1/2 «"+‘ 

= A, +1/2 h" + XC-t 1/2 A,_ 1/2 m” + 1 + XC^ 3/2 A i+3/2 u"+ 1 (2.21a) 
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He re C±_ 1/2 = C±(a»+/ /2 ). Next we take the absolute value of (2.21a) and 
use (2.18b) and the triangle inequality to obtain 

i + *C7 + i /2 + xc+ +1 /2 ]|a j+i/2 «“+ 1 | 

< |A )+1/2 u"| + XCt_ 1/2 |A ) _ 1/2 «"+ 1 | + XC7 +3/2 |A j+3/2 «"+‘| (2.21b) 
Rearranging terms, we get 

|A 3+1/2 «”+ 1 | < |A 3+l/2 «”| 

+ x C“)_ 3 / 2 |A : , +s / 2 ti"+ 1 | — C7+1/2 I Aj 4 .i/ 2 m" + i | 

- * Cf +i/ 2 |A 3+1/2 «"+ 1 | - C J + 1/2 |A 3 _ 1/2 t i "+ 1 | (2.21c) 

That is 

|Aj-(-i/2« n+1 | < |Aj+i/2« n | + MSj+i — Sy) (2.21d) 

where 

Si = C 3 -^ 1/2 |A 3+1/2 a’ , + 1 | - C+_ 1/2 |A 3 _ 1/2 «"+ 1 | (2.21e) 

Summing (2.21d) from j = — 00 to j = +00, we obtain (2.10); thus proving 
that our backward Euler implicit scheme is unconditionally TVD. 

We note that the TVD property does not depend on the particular form 
of the flux, but in general, is subject to a CF L like restriction (2.16). 

§2.4 Conversion to Second-Order Accurate Scheme 

Next, we want to briefly review the design principle behind the construction 
of second-order accurate TVD schemes. This is a rather general technique to 
convert a three-point first-order accurate (in space) TVD scheme (2.12) into 
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a five-point second-order accurate (in both time and space, or just space) 
TVD scheme of the same generic form. The design of high-resolution TVD 
schemes rests on the fact that the exact solution to (2.1) is TVD due to the 
phenomenon of propagation along characteristics, and is independent of the 
particular form of the flux f(u) in (2.1). Similarly, the first-order accurate 
scheme is TVD subject only to the CFL like restriction (2.16), independent 
of the particular form of the flux. Thus to achieve second-order accuracy 
while retaining the TVD property, we use the original TVD scheme with an 
appropriately modified flux (/ + g ), i.e., 

a? +1 + M/’+i/a - O = (2.22a) 



+ fj + 1 ~f 9] + ffj'+i 


Q( a j+l/2 H~ r fj+l/2)Aj+l/2u] 

(2.22b) 


where 


* 

lj + l/2 = ’ 


(fl'j+l — 9j)/&J + 1/2U 

0 


Aj-f 1/2 w 7 ^ o 

Aj-f 1/2U = 0 . 


(2.22c) 


The requirements on g are: (i) The function g should have a bounded 7 
in (2.22c) so that (2.22a) is TVD with respect to the modified flux (/ + g). 
(ii) The modified scheme should be second-order accurate (except at points 
of extrema; i.e., Uj = 1). In [2,3], Harten devised a recipe for g that 

satisfies the above two requirements. We will use this particular form of g for 
the discussion here. It can be written as 


gj = S • max [O, min(ay-j_i/ 2 | Aj-j-i/ 2 u|, S ■ (jj— 1/2 A j—i/2u)\ (2.22d) 

S = sign(Ay+i/ 2 u) 

with (Tj-f. 1/2 = a(aj-f. 1/2) and we choose 

a(z ) = - Q(z ) > 0 (2.22e) 

2 

for steady-state applications. It has the property that the steady-state solu- 
tion is independent of A£. Or, we choose 
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(2.22f) 


<7(2) = \( Q ( Z ) + x * 2 ) > 0 

for time accurate calculations. Note that if a(z) — ^{Q{z) + Xz 2 ), then 
(2.22) is second-order accurate in both time and space; see [3]. For transient 
calculations, second-order accurate in time is preferred. 

The form of g in (2.22d) satisfies the following relations [3]: 

(2.23a) 
(2.23b) 

(2.23c) 

Relation (2.23a) shows that the modified numerical flux (2.22b) is consistent 
with f(u). Relation (2.23b) shows that the mean-value characteristic speed 
lj+ 1/2 (2.22c) induced by the flux g is uniformly bounded. Relation (2.23c) 
implies that (2.22b) is second-order accurate in space. The form of g appears 
more complicated than it really is. The various test functions in (2.22d) can 
be viewed as an automatic way of controlling the numerical flux function so 
that (2.22) is TVD. 

The scheme (2.22) can be rewritten in the form (2.20) as 

«" +l -\ C -(a + 7 )”+, 1 /2 A J+1/2 u"+ l +XC+(a + 7)"l 1 1 /2 A ; ,'_ 1/2! /"+ 1 = u? 

(2.24) 

where C^a -f- ^±1/2 = ^j±i/ 2 )> i- e -> is now a function 

of (a -f- 7) instead of 0. The modified scheme (2.22) is of the same generic 
form as the original first-order scheme (2.17). Therefore (2.22) is an upstream 
differencing scheme with respect to the characteristic field (a + 7). Moreover, 
we have the following relation 

sign(a + 7) = sign(a) (2.25) 

for \z\ > 6, with z — a or (a 7). Hence (2.24) is also an upstream 
differencing scheme with respect to the original characteristic field a(u). 


9j = g{uj-i, Uj,u j+ i), g{u,u,u) = 0 
iTj-f-i/ 2 ! = IjTh-i 9j\/\ u j+i u j\ — ®i a j- 1-1/2) 


g = Axa(a)— + 0((Az) ) 
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Because of (2.23a), the numerical flux (2.22b) of the second-order accurate 
TVD scheme depends on four points; i.e., fj+ 1/2 — f( u j— i> u j> u j+i> u j+ 2 ), 
and thus (2.22) is formally a five point scheme. We note, however, that 

f{v,u,u,w) = f{u) (2.26) 

for all v and w. Hence, for practical purposes, such as numerical boundary 
conditions, (2.22) can be regarded as essentially a three-point scheme. 

We turn now to examine the behavior of TVD schemes around points of 
extrema, by considering their application to data where 

XI j— 1 ^ tij ■ — ^ tlj-^-2 (2.27) 

In this case gj = gj+i — 0 in (2.22d), and thus the numerical flux (2.22b) 
becomes identical to that of the original first-order accurate scheme (2.17b); 

A 

consequently, the truncation error of (2.22) deteriorates to 0((Ax) ) at j and 
; + l. This behavior is common to all TVD schemes. We conclude that, for 
a second-order accurate scheme to be TVD, it has to have a mechanism that 
switches itself into a first-order accurate TVD scheme at points of extrema. 
Because of the above property, second-order accurate TVD schemes are 
genuinely nonlinear, i.e., they are nonlinear even in the constant coefficient 
case. 

Extension of the one-parameter family of three- point TVD schemes (2.12) 
to second-order TVD schemes follows the same procedure except (2.22f) 
becomes 

°{z) = \q(z) - Mv - 1/2)* 2 (2.28) 


§2.5 Enhancement of Resolution by Artificial Compression 

Our technique to convert the first-order accurate TVD scheme (2.12) into 
a second-order accurate one is closely related to the concept of artificial 
compression (see [16,17]). 

Truncation error analysis shows that the first-order accurate scheme (2.12) 
is a second-order acc urate approximation to solutions of the modified equation 
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(2.29) 


u, + U = 


where a (a) is defined in (2.22e) or (2.28) 

We note that the CFL restriction (2.16) implies that a(a) > 0; thus the 
right-hand side of (2.29) is a viscosity term. Hence the first-order accurate 
TVD scheme (2.12) is a better approximation to the viscous equation (2.29) 
than it is to the original conservation, law. 

We obtain a second-order approximation to du/dt -j- df / dx = 0 by ap- 
plying the first-order scheme (2.12) to the modified flux (/ -f- g), where g is 
an approximation to the right-hand side of (2.29); i.e., 

9 = A2 ^( CT ^^) + °(( Ax ) 2 ) (2.30) 

The application of the first-order scheme to (/ -f- g) has the effect of 

cancelling the error due to the numerical viscosity to 0((Az) 2 ); thus g is 
an “anti-diffusion” flux. 

As to be expected, when we apply the first-order TVD scheme to (/ -f- 
(1 -j- u)g ), oj > 0, rather than to (/ -f- g), we find that the resolution of 
discontinuities improves with increasing u. This observation allows us to use 
the notion of artificial compression to enhance the resolution of discontinuities 
computed by the second-order accurate TVD scheme (2.22). This is done by 

increasing the size of g in (2.22d) by adding a term that is 0((Az) 2 ) in regions 
of smoothness, e.g., 


£/= (1 + u9j)gj, u > 0 

with 


) = jAj+i/ai/ - 
3 |^J + l/2 tt l + l/2U| 


(2.31a) 


(2.31b) 


Using g - (2.31) instead of g 3 makes the numerical characteristic speed more 
convergent, and therefore improves the resolution of computed shocks. Since 
$ = 0( As), this change does not adversely affect the order of accuracy of 
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the scheme. See [2] for more details. From numerical experiments, oj = 2 
seems is a good choice. 

We remark that applying too much artificial compression in a region of ex- 
pansion (i.e., divergence of the characteristic field a = df /du) may result in 
violation of the entropy condition. Hence when applying artificial compres- 
sion, one has to either turn it off in regions of expansion, or to limit the size 
of u) in (2.31a), say, by the value that makes (2.22) with (2.31a) third-order 
accurate (in regions of monotonicity). 


§2.6 Linearized Version of the Implicit TVD Scheme 

In order to solve for u n + l for this first- or second-order implicit scheme, 
we have to solve a set of nonlinear algebraic equations. The iterations needed 
to solve this system of nonlinear equations, say by Newton’s method at each 
time-step, may not compensate for the savings gained by the ability to take 
larger time-steps. 

To overcome this obstacle, we will present two ways of linearizing the 
implicit TVD scheme. The first method will preseve the conservative form 
of the differencing scheme but the resulting scheme may no longer be TVD 
or unconditionally TVD. The second method will destroy the conservative 
property but preserve its unconditionally TVD property. We will refer to 
the first method as the linearized conservative implicit (LCI) form, and the 
latter the linearized nonconservative implicit (LNI) form. The LNI is mainly 
useful for steady-state calculations, since the scheme is only conservative 
after the solution reaches steady-state. On the other hand, we have the 
advantage of stability of an unlimited CFL number. Note that the procedure 
of obtaining the LCI and LNI forms are applicable to both the first and second- 
order accurate Implicit TVD schemes. We will discuss the LCI and LNI for 
the second-order accurate one. To get the LCI and LNI for the first-order 
accurate TVD scheme, one simply sets q = 7 = 0 in the second-order form. 
Rewrite equation (2.22b) as 




4j + 1/2 W / 


4+l/2 U — Q{ a j+ 1/2 + 7j+1/2)4h-1/2 w 


n+ 1 


(2.32) 
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and use a local Taylor expansion about u n 

/ n + 1 - f n = a n (u n+1 - u n ) + 0(A£ 2 ) (2.33) 

where a = df/du. Applying the first-order approximation of (2.33) and 
locally linearizing the coefficients of A^x/gw in the second and third terms 
on the right-hand-side by dropping the time index from {n -f 1) to n, we get 
the LCI form 

“” + ‘ - |{ 5 + Ow - «“"+./ 2 + ^ + 1 / 2 ))^+ ,/ 2«" +1 

((*7—1/2 Q( a j— i j‘i + 7 7— — = u” (2.34a) 

where 



Ei = | -aj_, + - Q((a + 7)J_ 1/2 ) (2.35b) 

£2 = 1 + 2 — Pj+1/2 + < 2(( a + 7)?+l/ 2 ) — ^7- 1/2 

~f Q((° + 7)y_ i/ 2 ) (2.35c) 

£ 3=1 <+1 + /JJ +1/2 - <3((a + 7 )?+i/ 2) (2.35d) 
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and f’j+i / 2 is (2,22b)-(2.22e) calculated at the time level n. At this writing, 
we do not have any analysis to show that the LCI form is still TVD. 

The LNI form is obtained simply by replacing the coefficients ((?=*: ) rt + 1 
with (C^) n in (2.24), i.e., 


u"+ l —\C (a -f- 7)” +J/ f 2 A J+ i/ 2 w n " , " i +^C'" , "(< J + 7)j—i/2^j—i/2^ n+1 — u j 

(2.36) 

Since C ± > 0, it follows from (2.21) that (2.36) is unconditionally TVD. 

In delta form notation, (2.36) can be rewritten as 


1 — \C j+l/2 &j+i/2 + 


{u n+1 - u n ) 


X 




f j+1/2— f j-1/2 


(2.37a) 


where the left-hand-side equals 

dj \C j_^\j2&j+i/2d’ ~f" ^^^—1/2^3 — 1/2^ (2.37b) 

with dj = uj + 1 — uj, A j+i/ 2 d ~ d j+ 1 — dj, and Cf^ x / 2 == 
C±(a -f- 7 )J + i /2 s C±(aJ +1/2 + 7" +1 / 2 ). Rearranging terms, we get 

Eidj—i -f- E 2 dj -f- E^dj+i ~ — X[/y +1 / 2 — fj—i/ 2 ] (2.38a) 


with 


Ei = -\C+(a + 7)”_ 1/2 (2.38b) 

E 2 — 1 -f X[C , '“(a -f 7)” +1 / 2 -f- C + (a -j- 7)”_ 1 / 2 ] (2.38c) 

E 3 = -\C-(a + 7)"+i/2 (2.38d) 


18 



Again, f*+i / 2 is (2.22b)-(2.22e) calculated at the time-level n. It follows from 
(2.22b) and (2.37a) that the steady-state solution of (2.37) is 

(i) consistent with the conservation form, 

(ii) a spatially second-order accurate approximation to the steady-state of 
the partial differential equation, 

(iii) independent of the time-step At used in the iterations. 

Moreover, the iteration matrix associated with (2.38) is a diagonally 
dominant, tridiagonal matrix. Note that this linearized construction is not 
trivial, since the second-order method is a five-point scheme. Normally the 
matrix associated with (2.38) could have been a block pentadiagonal matrix. 
As mentioned before, (2.36) or (2.38) is not in conservation form and there- 
fore should not be used to approximate time-dependent solutions (transient 
solutions). However, it is a suitable scheme for the calculation of steady-state 
solutions. The rest of the paper will be devoted only to the LNI form. 


§3. Generalization to One-Dimensional Hyperbolic System of Conservation 
Laws 

At the present development, the concept of TVD schemes, like monotone 
schemes is only defined for nonlinear scalar conservation laws or constant 
coefficient hyperbolic systems. The main difficulty stems from the fact that, 
unlike the scalar case, the total variation in x of the solution to the system of 
nonlinear conservation laws is not necessarily a monotonic decreasing func- 
tion of time. The solution may actually increase at moments of interaction 
between waves. Not knowing a diminishing functional that bounds the total- 
variation in x in the system case, makes it impossible to fully extend the 
theory of the scalar case to the system case. What we can do at the moment 
is to extend the new scalar TVD scheme to system cases so that the resulting 
scheme is TVD for the u locally frozen ” constant coefficient system. To ac- 
complish this, we define at each point a u local" system of characteristic fields. 
This extension technique is a somewhat generalized version of the procedure 
suggested by Roe in [18]. 

Now, we briefly describe the above approach of extending the second-order 
accurate TVD schemes to hyperbolic systems of conservation laws 
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(3.1) 


SU . 9F(U) 
dt dx 


A{U) = 


9F 

dU 


Here U and F(U) are column vectors of m components and A(U ) is the 
Jacobian matrix. The assumption that (3.1) is hyperbolic implies that A(U ) 
has real eigenvalues a l (U) and a complete set of right-eigenvectors R l {U), 
l = 1, ..., m. Hence the matrix 


R(U) = (R'(U),..,R"‘(U)) (3.2a) 

is invertible. The rows L l (U), ...,L m (U) of R(U)~ X constitute an orthornor- 
mal set of left eigenvectors of A(U ); thus 

R- 1 AR = diag(a l ) (3.2b) 

Here diag(a*) denotes a diagonal matrix with diagonal elements a 1 . 

We define characteristic variables W with respect to the state U by 

W = R- X U (3.3) 


In the constant coefficient case (3.1) decouples into m scalar equations for 
the characteristic variables 


dw l 

~dt 


+ a< 


dw l 

dx 


= 0 , 


a 1 = constant 


(3.4) 


This offers a natural way of extending a scalar scheme to a constant coefficient 
system by applying it “scalarly” to each of the m scalar characteristic equa- 
tions (3.4). 

/ 

Let Uj+ 1/2 denote some symmetric average of Uj and Uj+\ (to be discussed 
later); i.e., 


U J+1/2 = Wj,Uj +t ) (3.5) 

Let Oy_|_ 1 / 2 , Rj-\- 1 /2> 1/2 denote the respective quantities of a 1 , R l , L l 

related to A{U 3 + 1/2). Let w l be the vector elements of W, and let a l J ^_ 1 / 2 = 

Wy_l_ 1 — ui l j be the component of Ay-j. 1/2*7 = Uj+\ — Uj in the /th charac- 
teristicjdirection; Le., 
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Aj-f 1/2*7 = ify+l/2<*j + l/2; «i+l/2 = Rj+l/2^j+l/2U (3.6) 

With the above notations, we can apply scheme (2.22) scalarly to each of 
the locally defined ( frozen coef ficient) characteristic variables of (3.1) as 
follows: 

Uj +1 + WV+,/2 - ^- 1 / 2 ) = Uj, (3.7a) 


Fj+ 1/2 = 2 ^ F i + ^V+i) 

1 m r 1 

+ Vi ^3 + 1 ~ Q( a j+ 1/2 + ^+l/2) Q! 5 + l/2 ^j+ 1/2 ( 3 -^b) 

z i=i L 

where 


5 • max 0, min(<rj +1/2 |a5. +1/2 |, 5 • ^_ 1/2 aJ_ 1/2 ) 
sign(aJ +1/2 ) 


, f (g l J+ 1 - *}■)/«$+ 1/ 2 «}+i /2 ^ 0 

7j+1/2 ~' 0 a* „ = 0 


(3.7c) 


(3.7d) 


Here <7 j_j_i/ 2 = cr i a j+ 1 / 2 ) where o{z) is (2.22e) and ot l J + l j 2 is (3.6). The 
corresponding in (2.31) for the added artificial compression term is 


= (1 + u‘ > 0 (3.7e) 

with 

e l __ I^J + l/2^ - Aj-i/ 2 m| ^ 

J |A i+1/2 u|-f IAj-j/smI 
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The cu l can be different from one characteristic field to another. 

Similarly, we generalize the linearized nonconservative implicit (LNI) form 
(2.37) to the system case by 


I —\J j+1 / 2 Aj+l/2 + 


(U n+1 - U n ) 


? 3+ 1/2 ~ Fj— 1/2 


(3.8a) 


or 


EiDj-i + E 2 Dj + E 3 D j+l = -X[F” +1/2 - f”_ 1/2 ] (3.8b) 

with 



Ei — 

(3.8c) 


F 2 = / + X[J- + 1 / 2 -h/t_ 1/2 ] 

(3.8d) 


Es — — X/y +1/ 1 2 

(3.8e) 

and 

jf +1/2 = fl" +1/2 diag(c±(a l + 7 ‘)" +i/ 2 )(R -')] +i/2 

(3.8f) 


D, = t/J+' - U" 

(3.8g) 

where the left hand side of (3.8a) is equal to 



Dj X/ J+1 y 2 Ayq-i/2D + X/^tlj^Ay—i^F 

(3.8h) 


with A y_{- i / 2 -F = Dy_|_i — Z>y. 


In the constant coefficient case where A(17) = constant, both (3.7) and 
(3.8) are TVD by construction. However, they are not identical; equation 
(3.7) is fully nonlinear while (3.8) is a version with a linearized left-hand-side. 

Note that the total variation for the vector mesh function U of the constant 
coefficient-case is defined as 
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(3.9) 


oo m 

TV{U) = E E K+1/2I 

j'= — co(= 1 


The particular form of averaging in (3.5) is essential if we require the scheme 
(3.7) for m — 1 to be identical to the scalar scheme of section 2, since we 
have to choose (3.5) so that Oj+ 1/2 * s same as the mean value in equation 
(2.17c). This can be accomplished by taking the eigenvalues a * + x / 2 and the 
eigenvectors f?y_|_ 1 y 2 in (3.2) to be those of A(Uj,'Uj+x), where A(Uj, Uj+i) 
is the mean value Jacobian. This matrix should satisfy 

(i) F{U) - F(V) = A(U, V){U — V) 

(ii) A(U, U ) = A(U) 

(iii) A(U, V ) has real eigenvalues and a complete set of eigenvectors. 

Roe in [18] constructs a mean value Jacobian for the Euler equations of 
gas dynamics of the form A{U,V) = Ai^U, V)), where ^{U, V) is some 
particular average. We will discuss Roe’s mean value Jacobian in the next 
two sections. 


§4. Applications to One-Dimensional Compressible Inviscid Equations of Gas 
Dynamics 

In this section we describe how to apply the implicit TVD scheme (3.8) 
to the compressible inviscid equations of gas dynamics (Euler equations). 
Included in this section are: (i) a detailed description of each of the terms 
of equations (3.5) to (3.8) for the Euler equations, (ii) Roe’s special form for 
A(Uj,Uj+ 1 ), (iii) an algorithm to compute U^ 1 , (iv) a description of the 
numerical example, and (v) a discussion of the numerical results. 


§4.1 Description of Algorithm 

In one spatial dimension, the Euler equations of gas dynamics can be 
written in the conservative form as 
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(4.1a) 


dU dF(U) 
dt + dx ~ 

where 


u = 

m 

, F = 

m 

™ 2 /p + p 


. e . 


_(e + p)m[p , „ 


(4.1b) 


Here U is the vector of conservative variables, F is the flux vector, and m = 
pu. The primitive variables are the density p, the velocity u, and the pressure 
p. The total energy per unit volume e, is defined as 


e — pe-\-pu 2 / 2 (4.1c) 

with e as the internal energy per unit mass. The pressure p for a perfect gas 
is defined as 


p = ( 7 — l)[e — m 2 /2p] (4. Id) 

where 7 is the ratio of specific heats and should not be confused with the 
13+1/2 in (2.22c) or l l j+i/2 in (3.7d). 

Let A denote the Jacobian matrix dF(U)fdU whose eigenvalues are 


(a l ,a 2 ,a 3 ) = (u — c,u,u-\- c) (4.2) 

where c is the local speed of sound. The right eigenvectors of A form the 
matrix R — (R l ,R 2 ,R 3 ) given by 


and 


R = 


1 1 1 
u — c u u-\- c 
L H — uc ^ u 2 H -f- uc, 

c 2 u 2 



i( 6 i + u f c ) 

i{-b 2 u - 1/c) 

£62 

R~ l = 

1 ~ h 

b 2 u 

— b 2 


M b 1 - u ! c ) 

h(~ b 2U + 1 /c) 

h b z. 


(4.3a) 


(4.3b) 


(4.3c) 
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with 



(4.3d) 

(4.3e) 


Let the grid spacing be denoted by Az such that x — j Ax. Using the same 
notation as in section 3 , the vector a of equation (3.6) is 


where 


a j + l/2 


' (aa — 66)/2 ' 

1/2 

= 

^■+i/ 2 P — aa 

, a ?+l/2. 


. (aa + 66)/2 . 


(4.4a) 


aa = 


7- 1 


c j + 1/2 L 


? ,2 

U j+ 1/2 

Aj+i/2e H — Aj-|-i/ 2P — Uj4-i/2Aj-j-i/2m 


66 = [A J+1 / 2 m — u J+1 /2&]+i/2p]/cj +1 / 2 


(4.4b) 

(4.4c) 


The simplest form of Uj+ 1/2 is 

U J+1/2 = (C/ i+1 + U i )/2 (4.5) 

Roe in [18] uses a special form of averaging that has the computational 
advantage of perfectly resolving stationary discontinuities. Roe’s averaging 
takes the following form 


u 


3 + 1/2 = = 


Du j+1 + Uj 


Hj+ 1/2 = — ^ 


D + 1 


ji 

S+1/2 

D 

H 


D + l 

(7- l)(i? J + l /2 - ltf+ l/2 ) 


— \J Pj + 1 / Pj 


+ I U 2 

• A 


(7 - l)/> 2 


(4.6a) 

(4.6b) 

(4.6c) 

(4.6d) 

(4.6e) 
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To use Roe’s averaging, we compute Uj+ i/ 2 , C7+1/2 ‘ n (4.2)-(4.4) by (4.6). 
In the numerical experiments for the one-dimensional test problem, we use 
Roe’s averaging. 

Given {/", for all j, we now list the operations needed to calculate 17-+ 1 
(assuming a fixed CFL number as input): 

(i) compute Uj = rrij/ pj and Pj 

(ii) compute Uj. 1-1/2 and Cj. j-1/2 from (4.5) or (4.6), calculate 

M = max(]uj-j-i/ 2 | + c j+i/ 2) 

3 

evaluate <4 '+i/ 2> / = 1, 2,3 by (4.4a), and define 

X = At/ Ax = fi/M 
where fi is the prescribed CFL number as input 

(iii) compute <4+1/2 b y ( 4 - 2 ); 4+1/2 b y ( 3Jc ) 

(iv) compute 7 j+i/ 2 b y (3.Td), and Fy+1/2 by (3.7b) and relation (4.3a) 

(v) compute C±(<4+ 1/2 + 7J+1/2) b y (2.18a) (with 0 = 4+ 1/2 + 7*+ 1/2 ) 

(vi) compute Jf^. l / 2 by (3.8f), and E\, E 2 and Es, by (3.8c)-(3.8e) 

(vii) solve the tridigonal system (3.8b) for Dj and then compute l/” 4-1 from 
(3.8g). 

§4.2 Numerical Example 

For the numerical experiments, a quasi-one-dimensional nozzle problem 
was selected. The governing equations for the nozzle problem can be written 
as 

f + «w,_» 

where 
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m/c 


U = 


' pK ' 


m/c 

, F = 

. e/c . 



( m 2 /p + p)k 



■ 0 ■ 

, H(U) = 

pfe 

J 

. 0 . 


(4.7b) 


with k, the area of the nozzle, a function of x. The nozzle we consider is a 
divergent nozzle [19] (figure 4. y with 


k(x) — 1.398 + 0.347tanh(0.8a; — 4) (4.7c) 

The steady flow condition were supersonic inflow, subsonic outflow with a 
shock. In all of the calculations the computational domain was 0 < x < 10. 
We used a very coarse mesh spacing of Ax = 0.5 (i.e., 20 spatial intervals), 
to evaluate the resolution of the scheme. 


Initial Conditions: 

We use linear interpolation between the exact steady-state boundary values 
as initial conditions. 


Analytical Boundary Conditions: 

We specified all three conservative variables p, u and e for the supersonic 
inflow, and the variable e for the subsonic outflow. 

Numerical Boundary Conditions: 

We used zeroth- or first-order space extrapolation to obtain the numerical 
boundary conditions for the unknown flow variables ( p and m) at the outflow 
boundary. Since the spatially second-order accurate TVD scheme is a five- 
point scheme, we also need the values of gj and 8j on both boundaries. For 
convenience, we will use zeroth-order space extrapolation for these values. 


§4.3 Discussion of Numerical Results 

All of the computations for the quasi-one-dimensional nozzle problem were 
done in single precision on the VAX 11/780 computer (a 6 digit machine). 
To illustrate the stability and/or accuracy of the LNI form of the implicit 
TVD scheme (3.8b), we compare in figure (4.2) the computed results with the 
explicit TVD scheme (forward Euler in time, obtained by setting Ei = 0, 
E 2 = and E& = 0 in equation (3.8b)), the first-order flux-vector splitting 
scheme [20], and a conventional implicit method using backward Euler in time 
and central spatial differencing with an added fourth-order explicit numerical 
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dissipation [21]. Figure (4.2a) shows the numerical result of the explicit TVD 
scheme. It took approximately 700 steps to converge to the steady-state with 
a fixed CFL = 0.8. Figure (4.2b) shows the computation of the LNI form 
of the second-order accurate implicit TVD scheme with CFL = 10 6 after 
25 steps. Both the explicit and implicit TVD schemes produce similar high 
resolution solutions (with 6 = 0.125 in equation (2.17d)). From numerical 
experiments, we found that artificial compression is not necessary for the 
one-dimensional problem. 

The convergence criterion is based on where all the terms on the right-hand 
side of equation (3.8b) are less than or equal to 10~ 4 . The implicit TVD 
method requires approximately triple the CPU time per time step more than 
the explicit TVD method but results in enhanced convergence rate. 

We have only used the nonconservative linearized delta form of the implicit 
TVD method for numerical experiments. A steady-state solution can be 
reached in 25-30 steps with CFL ranges from 10 6 — 10 7 . The steady-state 
solution profiles are independent of the CFL number. The number of steps for 
convergence monotonically decreases as the CFL number increases. However, 
the reduction in the number of steps is less pronounced for CFL number 
in the range from 10 3 to 10 6 . There are four primary factors affecting 
the convergence rate for CFL numbers higher than 10 3 : (1) the initial 

condition ( or the initial guess), (2) the numerical boundary conditions, (3) 
the interaction of nonlinear waves and (4) the machine accuracy of the VAX 
11/780. The influence of the above four factors will be the subject of a future 
investigation. 

Figure (4.2c) shows the converged solution by a conventional implicit 
method with CFL = 10. The oscillation near the shock is typical of a 
three-point central spatial difference scheme. The experimentally determined 
maximum CFL number is around 10 with the above initial and numerical 
boundary conditions. (We can improve the stability by adding an implicit 
second-order numerical dissipation term.) Figure (4. 2d) shows the converged 
solution by the first-order flux- vector splitting method [20] with CFL = 10 6 . 
Again, the steady-state is reached after 25 — 30 time steps. The solution is 
very smeared but is independent of the CFL number. From the results, we 
can see a definite improvement in shock resolution by the TVD schemes over 
the conventional methods. The implicit TVD scheme requires approximately 
80% more CPU time per time step than the conventional implicit method. 
Moreover, there is a dramatic increase in convergence rate of the implicit 
TVD scheme over the explicit TVD scheme. 
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§5. Applications to Two-Dimensional Compressible Inviscid Equations of Gas 
Dynamics 

In this section we describe how to formally extend the one-dimensional im- 
plicit TVD scheme to an Alternating Direction Implicit {ADI) version for 
the two-dimensional compressible inviscid equations of gas dynamics (Euler 
equations). This is a formal extension of the TVD scheme from one dimen- 
sion to two dimensions. At the present state of development, there is not 
yet a similar theoretical analysis of the TVD properties for two-dimensional 
hyperbolic equations. Included in this section are: (i) a detailed descrip- 
tion of the analogue of each of the terms of equations (4.2) to (4.6) for the 
two-dimensional Euler equations, (ii) a discussion on the extension of the ex- 
plicit and implicit TVD schemes to two-dimensions, (iii) Roe’s form for the 
Jacobian matrices A{Uj t k, and B{Uj t k, E/y+i,*), (iv) an algorithm to 

compute Ufy 1 by the ADI approach, (v) a description of a numerical ex- 
ample, and (vi) a discussion of the numerical results. 

§5.1 Numerical Fluxes in Two-Dimensions 

In two spatial dimensions, the Euler equations of gas dynamics can be 
written in the conservative form as 


where 


dU dF{U) dG{U ) 

dt dx dy 


(5.1a) 


U == 


-p- 


m 


n 

m 


m 2 /p + p 


nu 

n 

, F = 

mv 

, G = 

n 2 /p + p 

- e - 


.(e + p)m/p. 


.(e + p)n/p. 


(5.1b) 


with m = pu and n = pv. The primitive variables are the density p, the 
velocity components u and v, and the pressure p. The total energy per unit 
volume e, is related to p by the equation of state for a perfect gas 


P = (7- 1) 


e 


(m 2 + n 2 )' 

' 2 ~p . 


(5.1c) 


where 7 is the ratio of specific heats and should not be confused with the 
7j+i/2 in (2.22c) or l l j+1/2 in ( 3 - 7d )- 
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Let A denote the Jacobian matrix dF(U)/dU whose eigenvalues are 


(<4, al,al,a*) = (u — c,u, u + c,u) 


(5.2) 


where c is the local speed of sound. The right eigenvectors of A form the 
matrix R x = {R\,R 2 ,R*,R\) given by 


R x = 


1 

u — c 

V 


1 

u 

V 


1 0 

u-\- c 0 

v 1 


\_H — uc ( u 2 -f- v 2 )/2 H -f uc vj 


(5.3a) 


where 


and 


H 


1- 1 


+ 


u 2 + tr 


R 


— i 


with 


’i(*l + «/c) i(— &2«— l/c) i(~ 6 2 v) ^2 


6 2 u 


6 2 V 


-6c 


4(*1 — u/c) i(— 6 2 u + 1/c) £(— 6 2 v) £6 2 


— v 


0 


0 J 


(5.3b) 


(5.3c) 


61 = 6 2 


(u 2 -f ^ 2 ) 


6 2 = 


7 - 1 


(5.3d) 

(5.3e) 


Let the grid spacing be denoted by Ax and Ay such that x = j Ax and 
y = kAy. Using the same notation as in section 4 , the vector a of equation 
(3.6) for the ^-direction (omitting the k index) is 
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a J+l/2 


(aa — 66)/2 

a H-i/2 


Ay+i/ 2 p — aa 

"f+1/2 


(aa -j- 66)/2 

AW 


-A‘+l/2 n ~ v j+l/2&j+l/2P. 


where 


(5.4a) 


aa — 


1- 1 
C H*l/2 


A'-f 1/26 -}- 


U H~l/2 + V ?+l/2 
+ 1/2 P 


— *fH-i/2A'-H/2 m — Vj4. 1 / 2 Ay +1 / 2 n 
bb = [A J+ i/ 2 m — «y + i/ 2 Ay + i/ 2 p]/cy +1 /2 


(5.4b) 

(5.4c) 


Similarly, let B denote the Jacobian matrix dG(U)JdU whose eigenvalues 
are 


(«J. <4> “5> a J) = (»-«. V, V + C,v) 


(5.5) 


where c is the local speed of sound. The right eigenvectors of B form the 
matrix R y = (i?J, R 2 , R z y , R*) given by 


and 


Ru — 


1 

u 

V — c 


1 

u 

V 


H — vc ( u 2 -j- v 2 )/2 


1 

u 

V -j- c 
H + vc 


O' 

1 

0 

u 


(5.6a) 



h(t>i + v/c) £(- b 2 u ) 
1 — 6i b 2 u 

i(*i ~ v / c ) h(~ b 2u) 

— u 1 


£(-6 2 v — 1/c). £& 2 ' 
b 2 v —b 2 

b 2 v 1 / c) £& 2 

0 0 . 


(5.6b) 
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with 61 and 62 defined in equations (5.3d) and (5.3e). The vector a for the 
y-direction (omitting the j index) is 


(cc — dd)/2 

A*-fl/2/> — cc 

(cc + dd)/2 

where 


(5.7b) 
(5.7c) 

As mentioned previously, the simplest form for ££,+ 1 / 2 ,* is 

££ 7 + 1 / 2 ,* = (Uj+i,k + Uj,k)/2 (5.8) 

Roe's special form of the averaging in the ^-direction is: 


W ft + l/2 + V * + l/2 


&k+l/2p 


• f 

~ WA; + l/2A/c-}_i/2m — Vk-\-l/2^k + l/2n 
dd = [Afc_|_i/ 2 rc — Vfc^_i/ 2 Afc^.i/ 2 p]/cfc-j.i /2 


(5.7a) 


Q l + l/2 

Q k+ 1/2 _ 

«S+i/a 

Q Hl/2. 
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ttj+1/2,* = 
Vj+l/2,k = 
Hj+l/2,k = 
4-H/2,* = 

D = 
H = 


DUj+l,k H~ u j,k 

D+l 

Dvj+I,k + Vj,k 

D+l 


DHj+ i,k "f H 3 ,k 


D+ 1 


(7- 1) 


Hj+l/2,k — ~(Uj +l / 2 k + tf+l/2, k) 


\/ 'Pj+iyk/ Pj, 


IP 


h-l)p + 2 ( “ 2+ ^ 


(5.9a) 

(5.9b) 

(5.9c) 

(5.9d) 

(5.9e) 

(5.9f) 


Therefore to use Roe’s averaging for the z-differencing, all one has to do is 
compute Uj+i/ 2 ,k, Uj+ 1 / 2 and c j+i/2,k in (5.2)-(5.4) by (5.9). Similarly, we 
can obtain the Roe’s averaging for Uj >k +i/2, v j,k+ 1/2, and Cj tk + 1 / 2 - In the 
numerical experiments for the two-dimensional test problem, we used Roe’s 
averaging. 

The two-dimensional form of numerical fluxes (3.7b) for the Euler equations 
of gas dynamics (5.1) can be written as 


~ n 


Fj+i/2,k = iinv i ,k)+nu i +i,k'A 


1 

+ o 2J ff 3 + gt J+l ~ Q( a j+l/2 + 7>+l/2) a i+l/2 


Rl J + 1/2 


(5.10a) 


~ n 


^>,* + 1/2 ~ «[ G ( U 3,k) + G (Uj,k- f-l)] 


1 4 r 

+ 9 S 0* + 9k + 1 — ^( fl *+l/2 + 7*-|-i/2) a *+l/2 


f=i L 


^i-(-i/2(^*10b) 


33 



with 9j-^~ i and 9k+i corresponding to (3.7c) evaluated at the appropriate 
“ locally frozen ” Ith characteristic speed and /th characteristic variable in 
the x- and y-direction, and Q(z) is defined in (2.17d). The otj+ 1/2 and a k+i /2 
are defined in (5.4) and (5.7). Also, the following notations have been used 


a j+ 1/2 = ( a x)j+l/2,k 

(5.10c) 

a k + l/2 ~ ( a y ) j > k + 1/2 

(5.10d) 

1/2 = 1/2, k 

(5.10e) 

Tk+l/2 — hy)j )k+ 1/2 

(5.10f) 


where a l x is defined in (5.2), a l y is defined in (5.5) and 


1 f (?5+i — ?5)/ Q i+i/2 “ 3 + 1/2 ^ 0 

7 ,+./ 2-| 0 


1/2 — 0 


(5.10g) 


lk+ 1/2 




[ (9 l k + 1 — 9k)/ a k + 1/2 a k + l/2 7^ 0 

0 Q i4-i/2 == 0j 


(5.10h) 


Here, it is understood that the scalar values and the vector R l in the 
summation in equations (5.10a) and (5.10b) are values of (5.2)-(5.9) evaluated 
at the coresponding x- and ^-coordinates. For simplicity, we omitted the k 
index inside the summation sign of equation (5.10a), and omitted the j index 
inside the summation sign of equation (5.10b). 


§5.2 Extension of the Explicit TVD Scheme by the Fractional Step Method 

In this subsection, we are going to review the implementation of the ex- 
plicit TVD scheme to two dimensions by the fractional step (time splitting) 
method. We will also give a discussion on the use of the artificial compression 
term. Later, we will show a comparision between the explicit and implicit 
TVD schemes. 
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The explicit TV D scheme can be implemented in two space dimensions by 
the method of fractional steps as follows: 


V)* = 0J,» - - f'i-Mv) = UVjj! 

i h * At ~ * ~ * * 

Vft 1 = U i,k ~ - G 3 ,*-i/ 2 ) = 

that is 

C/"+ l = LyUUl, (5.12) 

where F* +1 j 2t k and * 4 . 1/2 are defined in (5.10). 

In order to retain the original time accuracy of the method, we use a Strang 
type of fractional step operators, namely 

c/"t 2 = UL y L y L x Ul k (5.13) 

or for steady- state calculations, we use 


(5.11a) 

(5.11b) 


w $ 2 = LlL,L,L,...L,Llu° Jtk (5.14) 

where L* denotes the operator with the time step equal to At/2. 

We may enhance resolution the same way as discussed in section 2.5 by 
increasing the value of 7 y-j-i /2 ( an d 7* 4 - 1 / 2 ) equation (5.10g) (and (5.10h)). 

To accomplish this, we increase the size of the corresponding < 7 * in (3.7c), for 
example, by multiplying the right-hand-side of (3.7c) by [1 -j- i.e. 

0, = [1 + }g‘j, u‘>0 (5.15a) 

with 


35 



(5.15b) 


< _ l Q J+l/2 ~ — 1/2 1 

^ \ Q j+l/2\ H" l Q j— 1/2I 

where is defined by (5.4a). Similarity, we can obtain <7* = [1— |-a ^0 l k ]g l k 

for the ^/-direction. 

A preliminary experiment in two-dimensional calculations using the explicit 
TVD scheme by the fractional step approach indicated a need for such an 
enhancement mechanism [2]. The artificial compression term is especially 
needed for the linearly degenerate characteristic field, i.e., a x — u and a y = 
v. We will use this form of the g function for our two-dimensional numerical 
experiments. 


§5.3 Extension of the Implicit Scheme by the Alternating Direction Implicit 
( ADI) Method 

The two-dimensional linearized nonconservative implicit (LNI) form of (3.8) 
for the Euler equations of gas dynamics (5.1) can be written as 


1 ~ J j + 1/2,* Ah- 1/2,* + — 1/2, fc 


-WK- k+I/2 A i>t+1/2 + \»K+ t — 1/2^7,* — 1/2 


(U n+1 - U n ) 


= -X 5 


~ n 


F j+l/2,k — F j-\/2,k 


- X* 


S,*+ 1/2 ~ U j,k- 1/2 


(5.16a) 


with X* = At/ Ax, \ y = At/ Ay, where 


Jf+i/2,k — (-Ra;diag(C^ : )i? ;z . . l j 2> k (5.16b) 

*&+./ 2 = fording {C?)R~% k+l/2 . (5.16c) 


and 
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(C± )y + i/ 2l * = 7j\Q( a x + l l x)±( al x + 7 l x)]"+i/2,k (5.16(1) 

( C ^ )i,fc+ 1/2 = + ^±K + ^C*+l/2 ( 5 - 16e ) 

/ = 1 , 2 , 3, 4 


It is well known that solving the two-dimensional implicit difference equa- 
tion (5.16) is very costly. This leads to the popularity of using the alternating 
direction implicit (ADI) method to solve gas dynamics problems. Formally, 
we can write an ADI form of (5.16) as 


I 


I 


J j+l/2,k^3 + l/2,k + l/2,lfcA?— 1/2, fc 


D* = 


- n 




F i+l/2,fc ~ F j-l/2,k 


-\ y 


Gj,k+ 1/2 Gj,k- 1/2 


\ y K j k+1/2 Aj- fc+i/2 + X s ' J ftr+ fc _ 1/2 A ;iA _i/2 


D = D* 


u n+l =U n +D 


(5.17a) 

(5.17b) 

(5.17c) 


Given for each O', A:), we now list the operations needed to calculate 
Ufy 1 by the ADI method (assuming a given CFL number): 

(i) compute u j>k — m J)k /pj tk , v j]k = n jtk /pj,k and p jik 

(ii) compute Uj+\/ 2 ,k, Vj+\/ 2 ,k and Cj. 1 - 1 / 2 , a from (5.8) or (5.9), calculate 


Afx — max(|Uy_f-i/2,&| “1“ Cj- {-1/2, fc) 

3 )* 

My = max(|^_ H/2> *| + Cj+i/ 2 , k ) 

3,K 

evaluate Oij +1 / 2 , l = 1, 2, 3, 4 by (5.4a), (remember that the k index is omitted 
from the equation), and define 
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\ x = At/ Ax = p/max (M z ,M y ) 

j,k 

where n is the prescribed CFL as input 

(iii) compute {a l x ) j+1/2k by (5.2), ff$ + i/ 2f * by (3.7c), (with (ai) i+1/2|fc as 

a j+l / 2 (3.7 c)) 

(iv) compute (7^.) J+1 / 2) fc by (5.10g), and Fj+ i/ 2) * by (5.10a) and relation 
(5.3a) 

(v) compute C± by (5.16d) (with 2 = {a l x ) j+l / 2 ,k+(7z)j+i/2,k in equation 
(2.18a)) 

(vi) compute Jf^. l / 2k by (5.16b) 

(vii) compute u jili+l / 2 , ^j,k+ 1/2 and Cj jk+x / 2/ calculate 

M x = max(!u ilfc+1 / 2 | + c j)k + 1/2) 

M y = max(|i>y ( *-pi/ 2 | -f- c Jr ,fc_|_i/ 2 ) 

evaluate 1/2^ — 1; 2 , 3, 4 by (5.7a) (remember that the j index is omitted 
from the equation) and define 


\ y = At/ Ay — /z/max (M z , M y ) 

3,k 

(ix) compute (a l y) j k+1/2 by (5.5), g l jjk+l/ 2 by (3.7c), (with {a l x ) jk+1/2 as 
4+ 1/2 in ( 3 - 7c )) 

(x) compute (7^. k+l j 2 by (5.10h), and G j>k+ x / 2 by (5.10b) and relation 
(5.6a) 

(xi) solve equation (5.17a) for D* 

(xii) compute C± by (5.16e) (with 2 = {a l y ) j k+1/2 + (^y) J)Jk+1 / 2 in ec l ua " 
tion (2.18a)) 
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(xiii) compute K^ k+l f 2 by (5.16c) 

(xiv) solve equation (5.17c) for Ufy 1 . 

We now turn to the discussion of the influence of the approximate fac- 
torization error on rate of convergence for the ADI method. It is well known 
that contrary to the fact that the convergence rate for the unfactored scheme 
is directly proportional to the CFL number, the convergent rate for a factored 
scheme has a definite maximum. That is, the iteration count grows rapidly 
when the calculation is carried out away from an optimal time- step. Recently, 
Abarbanel et. al. [22,23] have made some detailed analyses of the influence 
of the approximate factorization error on the efficiency of convergence rate as 
a function of CFL number for a conventional ADI scheme. We feel that, in 
our case, this phenomenon is related to the fact that the ADI version ceased 
to be TVD (in the scalar case) for large enough CFL number. At this time, 
however, we do not have a rigorous analysis for the dependence of the rate 
of convergence on the CFL number. Here, we extend (5.16) to an ADI form 
for the purpose of numerical experiments. 


§5.4 Numerical Example 

In order to examine the applicability of the new method for two- 
dimensional shock calculations, we consider a simple inviscid flow field 
developed by a shock wave reflecting from a rigid surface (figure (5.1)). The 
steady-state solution can be calculated exactly and thus can aid us in evaluat- 
ing the quality of the numerical method. Figure (5.1) shows the indexing of 
the computational mesh. 

Initial Conditions : 

Initially, the entire flow field is set equal to the freestream supersonic inflow 
values plus the analytical boundary conditions as described below. 

Analytical Boundary Conditions : 

The boundary conditions are given as follows: (a) supersonic inflow at 
j = l f lc = 1 which allows the values to be fixed at freestream 

conditions; (b) prescribed fixed values of Uj t k at k — K, j = 1, ..., J which 
produce the desired shock strength and shock angle; (c) supersonic outflow 
at j = J, k = 1, ..., K] (d) a rigid flat surface at k = 1, j = 1, ..., J which 
can be shown to be properly represented by the condition v = 0, with the 
additional condition dp/dy — 0 at k .= 1 from the normal y-momentum 
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equation. 




Numerical Boundary Conditions : 

The supersonic outflow values Uj ik , k = 1 — 1 are obtained by 

zeroth-order extrapolation, i.e., 


U J)k = Uj- 1)k , k = - 1 (5.18) 

The values of p Jt i, m^i on the rigid surface with j = 1, ..., / are also obtained 
by zeroth-order extrapolation, i.e., 


Pj.i — Pj, 2 

j = (5.19) 

™J,1 = W i,2 

Three different methods were used in approximating dp/dy = 0 to get 
e j , l > J = 1 > ’ * • > J 

(a) normal derivative of e equal to zero (first-order) 

fijji — ej } 2 (5.20a) 


(b) normal derivative of e equal to zero (second-order) 


e j, 1 = ( 4e j,2 ~ Cj, 3)/3 


(5.20b) 


(c) normal derivative of p equal to zero (second-order) 


Pj.i = (%,2 - Pj, s)/3 


(5.20c) 


e i,i 


Pj.i 


+ 


m 




(7-1) 2/7, ! 


(5.20d) 


Equation (5.20a) together with equation (5.19) is an approximation to 
Pj t i = Pj' 2 ', i e., a first-order approximation for the normal derivative of 
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p equal to zero. Equation (5.20b) together with equation (5.19) is an ap- 
proximation to equation (5.20c). We use equations (5.20a) or (5.20b) for the 
implicit method mainly because of their ease of application with implicit 
numerical boundary conditions. From the numerical experiments, we found 
that equation (5.20b) and (5.20c) produce better numerical solutions near the 
wall than equation (5.20a). 

Since this implicit TVD scheme is a 5-point scheme (in each spatial direc- 
tion), we also need the values of g likt g Jtk , g j} i, g jt K, 0i,&, 0j tk , 9 3>1 , 6 3>K , for 
j = 1, ..., J or k = 1, ...,K . For convenience, we will set 

9i,k = 92, k 0i,* = 02,* 

gj,k — gj— i,* 0 /,* = 0 /— i,* 

9 j,l = 9 j,2 0J,1 = 0J, 2 ( 5 - 21 ) 

9j,K ~ 9j,K - 1 0J.K = 6j,K - 1 


§5.5 Discussion of Numerical Results 

The purpose of these numerical experiments is three fold: 

(i) To test the performance of the ADI form of the implicit TVD scheme 
(5.17) on the shock reflection problem. For reference purposes, we denote 
scheme (5.17) as STVD. 

(ii) To test the performance of scheme (5.17) with qj. == 7^ = 0 in (5.16b) 

and (5.16c) (i.e., first-order spatial difference for the implicit operator) on the 
shock reflection problem. For reference purposes, we denote this method as 
FTVD. 

(iii) To test the effect of the choice of a and Q function in (2.22), and the 
artificial compression term in (5.15) on the convergence rate and resolution 
of the above two-dimensional problem. 

0 

In all of the numerical experiments, the incident shock angle ip was 29 
and the freestream Mach number Moo was 2.9. The computational domain 
was 0 < x < 4.1, and 0 < y < 1, with a uniform grid size of 61x21. The 
appropriate analytical boundary conditions were applied along the boundaries 
of the domain. Only pressure contours and pressure coefficients will be 
illustrated. Here, the pressure coefficient is defined as c p = ( p/Poo — 1) 
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with 7 the ratio of specific heats, and the freestream pressure. The exact 

o 

minimum pressure corresponding to ip = 29 and M oo = 2.9 is 0.714286 and 
the exact maximum pressure is 2.93398. The exact pressure solution and 
the computation domain is shown in figure (5.2). Forty-one pressure contour 
levels between the values of 0 and 4 with uniform increment 0.1 were used for 
the contour plots. The pressure coefficient was evaluated at y = 0.5 for 0 < 
x < 4.1. All of the computations for the two-dimensional shock reflection 
problem were done in single precision on a CDC 7600 computer. 

Comparison of Method 

Pressure contours and the pressure coefficients evaluated at y = 0.5 are 
shown in figure (5.3) for four different methods. Figure (5.3a) shows the 
numerical result of the explicit TVD scheme by the fractional step method 
(5.14) with 


Q( z ) = z 2 -J- 1/4 (5.22a) 

0 ( 2 ) =1/8 (5.22b) 


It took approximately 350 steps to converge with a fixed CFL = 0.8. The 
average smearing of the shocks is two points. Figure (5.3b) shows the steady- 
state solution of the conventional implicit method [21] with CFL = 1.0 after 
approximately 600 steps. The oscillations are spread over 8 grid points. The 
experimental maximum CFL number for this conventional implicit method 
is around 1.5. Figure (5.3c) shows the numerical result of the ADI form of 
the implicit TVD scheme (5.17) ( STVD method with 6 = 0.125 and <7 in 
equation (2.22e) ) after 300 steps with CFL = 3. The experimental maximum 
CFL number for this method is 5 under the CFL sampling sequence of 
(1,2,3,5,10,20). Figure (5.3d) shows the numerical result of the FTVD 
method (with 6 = 0.125 and o in equation (2.22e) ) after 60 steps with 
CFL = 6. The average smearing of the shocks is 2 points. The convergence 
rate of the FTVD method is far better than the STVD method. It was found 
that the optimal CFL (with A t fixed) for fastest convergence rate with the 
FTVD method is between 5 and 10. 

In order to show that the convergence rate is not a monotone decreasing 
function of CFL number for the FTVD method, figure (5.4a) illustrates the 
same method as figure (5.3d) after 250 steps with CFL = 50. The solution 
is not quite convergent yet. Figure (5.4b) shows the same method after 250 
steps with CFL — 500 (We ran for another 400 steps, without reaching 
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steady-state). Figure (5.4c) shows the same method after 250 steps with CFL 
= 1000 . 

A primary factor affecting the stability and convergence rate of the ADI 
form of the implicit TVD scheme is the approximate factorization error. This 
will be a subject of further investigation. 

The Choice of a and Q Functions, and the Artificial Compression Term 

The primary difference in shock resolution between the explicit and implicit 
TVD methods of figures(5.3a) and (5.3c) is the result of using a different a 
and Q. Equations (2.22e) and (2.17d) are used for the implicit scheme and 
equation (5.22) is used for the explicit scheme. We have found that the 
shock resolution is somewhat degraded by the explicit TVD scheme if a 
in (2.22e) is used instead of (5.22). The particular choice of a in (2.22e) is 
especially appropriate for the ADI method since the steady- state solutions 
are independent of At. 

All of the TVD schemes were operators with an artificial compression term 
(5.15a) with oj 1 — 2 ,1 — 1,2, 3, 4 except for the FTVD method. For the 
FTVD method, we used w 1 = w 3 = 1 and u 2 = w 4 = 2. We found that 
with this set of u l , we get a faster convergence rate. 

Approximate CPU Time and Actual Implementation of the Numerical 
Flux Functions 

The conventional implicit method requires approximately twice the CPU 
time per time step as the explicit TVD scheme. The implicit TVD method 
requires approximately 2.5 times more CPU time per time step than the 
conventional implicit method. 

In actual implementation of the explicit and implicit TVD methods into a 
computer code, the following form of the numerical flux F i+i* k (similarily 

for fc+ 1 / 2 ) was used instead of (5.10) 


F j+ U2,k = ^nu jt k)^F(U j+l , k )) 

+ 0 Sf^'+l/2 a i'-i-l/2(^ H" 0J + 1 ) — Q( a \+ 1/2 + 7j'+l/ 2 ) a i+l/2 R j+ 1/2 

(5.23a) 


'<=i' 
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with 

£}•+!/: 2 = 1 + u l ^ ^y+i) (5.23b) 

where is defined in equation (5.15b), and 


= 5 • max 


0,min(|a' + 1/2 |,S-a‘_ 1/2 ) 


5 = sign(a‘ +1/2 ) 


and 


(5.23c) 


'lj+l/2 = 


cl n l 

<3 + 1/2° j+l/2 


i (Sj + 1 ~ ^)/ a i+i/2 
0 


a j+l/2 7^ 0 

a J+l/2 = O' 


(5.23d) 


That is, the o l J+1/2 has been taken out from the definition of g< in equation 

(3.7c) for simplicity. Furthermore, the artificial compression is incorporated 
into the definition of the numerical flux function. 


§6. Concluding Remarks 

The nonlinear, spatially second-order accurate, unconditionally stable im- 
plicit TVD scheme in a “linearized nonconservative” form has been applied 
to obtain steady-state solutions for the one-dimensional compressible invis- 
cid equations of gas dynamics. This linearized form of the implicit TVD 
scheme is only conservative after the solution reaches steady-state. Numerical 
experiments for a quasi-one-dimensional nozzle problem show that the ex- 
perimentally determined stability limit correlates exactly with the theoretical 
stability limit for the nonlinear scalar hyperbolic conservation laws. Steady- 
state solution can be reached in approximately 25 steps. 

We have formally extended the second-order accurate implicit TVD scheme 
by an ADI method for two-dimensional calculations. Numerical experiments 
with the ADI form (5.17) for a shock reflection problem show the gain in 
efficiency is not as pronounced as the one-dimensional counterpart. 
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A spatially first-order accurate left-hand side (i.e., by setting 7^. = 7*, = 0 
in (5.11b) and (5.11c)) of the ADI form provides better stability and a faster 
convergence rate. A steady-state solution can be reached in approximately 60 
steps for a two-dimensional shock reflection problem. Numerical experiments 
also show that the rate of convergence is very sensitive to the CFL number. 
The iteration count grows rapidly when the calculation is carried out away 
from an optimal time- step. 

More rigorous analysis of the influence of approximate factorization error 
on the stability and efficiency of the method is needed. This will be the 
subject of a future investigation. 
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(a) Explicit TFD method. 


(b) Implicit TVD method. 



(c) Conventional implicit method (^) First-order implicit flux- vector 

splitting method. 


Fig. 4.2 Density distribution: supersonic inflow, subsonic outflow; 
20 spatial intervals. 


50 






y 



Fig. 5.1 Indexing of computational mesh 
for shock reflection problem. 
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Fig. 5.2 The exact pressure solution for 
the shock reflection problem. 
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PRESSURE CONTOURS 


PRESSURE COEFFICIENTS 



(a) Explicit TVD method, 350 steps, CFL — 0.8. 



(b) Conventional implicit method, 600 steps, CFL = 1. 



(c) Implicit TVD method ( STVD ), 300 steps, CFL = 3. 



(d) Implicit TVD method (FTVD), 60 steps, CFL = 6. 

Fig. 5.3 Pressure contmirs and pressure coefficients-for the shock 

reflection problem. 
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(b) 250 steps, CFL = 500. 


y 



(c) 250 steps, CFL — 1000. 


Fig. 5.4 Pressure Contours for three different CFL numbers. 
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